"""
    BuddyCheck

Spatial consistency quality control through buddy checking. This module identifies
isolated outliers by comparing each observation with nearby observations of the
same type.

# Mathematical Foundation

## Buddy Mean
For an observation y_i with buddies B_i:
```
ȳ_B = (1/|B_i|) Σ_{j∈B_i} y_j
```

## Distance-Weighted Buddy Mean
```
ȳ_B = Σ_{j∈B_i} w_j * y_j / Σ_{j∈B_i} w_j
where w_j = exp(-d_ij / L)
```

## Buddy Check Statistic
```
d_buddy = (y_i - ȳ_B) / √(σ_i² + σ_B²/N)
```

where:
- σ_i: observation error
- σ_B: buddy standard deviation
- N: number of buddies

## Decision Rule
Reject if |d_buddy| > threshold (typically 2-3σ)

# Buddy Search Parameters

**Horizontal radius**: 50-300 km (variable type dependent)
**Vertical radius**: 50-500 hPa (for upper-air obs)
**Temporal window**: 30-180 minutes
**Minimum buddies**: 3-5 observations

# Special Cases

**Isolated stations**: Pass if insufficient buddies
**Coastal regions**: Use land/sea mask
**Topography**: Adjust for elevation differences
**Different platforms**: Platform-specific weighting
"""
module BuddyCheck

using LinearAlgebra
using Statistics
using Dates
using ...ObservationTypes
using ...GSICoreAnalysis: AbstractStateVector

# Export types and functions
export BuddyCheckConfig, BuddyCheckResult
export perform_buddy_check, find_buddy_observations, compute_buddy_statistic
export distance_weighted_mean, apply_elevation_correction

"""
    BuddyCheckConfig{T<:AbstractFloat}

Configuration for buddy check quality control.

# Fields
- `horizontal_radius::T`: Horizontal search radius in meters
- `vertical_radius::T`: Vertical search radius in Pascals
- `temporal_window::T`: Temporal window in seconds
- `min_buddies::Int`: Minimum number of buddy observations required
- `rejection_threshold::T`: Threshold for buddy check rejection (sigma)
- `distance_weighting::Bool`: Use distance-based weighting
- `elevation_correction::Bool`: Apply elevation correction for surface obs
- `platform_weighting::Bool`: Weight by platform quality
- `variable_specific_radius::Dict{ObservationVariable,T}`: Custom radius per variable
"""
struct BuddyCheckConfig{T<:AbstractFloat}
    horizontal_radius::T
    vertical_radius::T
    temporal_window::T
    min_buddies::Int
    rejection_threshold::T
    distance_weighting::Bool
    elevation_correction::Bool
    platform_weighting::Bool
    variable_specific_radius::Dict{ObservationVariable,T}

    function BuddyCheckConfig{T}(;
        horizontal_radius = T(100_000),    # 100 km default
        vertical_radius = T(5000),          # 50 hPa default
        temporal_window = T(3600),          # 1 hour default
        min_buddies = 3,
        rejection_threshold = T(2.5),
        distance_weighting = true,
        elevation_correction = true,
        platform_weighting = true
    ) where T

        # Variable-specific search radii
        var_radius = Dict{ObservationVariable,T}(
            # Conventional observations - tighter radius
            TEMPERATURE => T(100_000),         # 100 km
            PRESSURE => T(150_000),            # 150 km
            SURFACE_PRESSURE => T(150_000),    # 150 km
            WIND_U => T(100_000),             # 100 km
            WIND_V => T(100_000),             # 100 km

            # Satellite observations - wider radius
            BRIGHTNESS_TEMPERATURE => T(200_000),  # 200 km

            # Specialized observations
            GPS_REFRACTIVITY => T(300_000),    # 300 km
            RADAR_REFLECTIVITY => T(50_000),   # 50 km (shorter for radar)
        )

        new{T}(
            horizontal_radius,
            vertical_radius,
            temporal_window,
            min_buddies,
            rejection_threshold,
            distance_weighting,
            elevation_correction,
            platform_weighting,
            var_radius
        )
    end
end

BuddyCheckConfig() = BuddyCheckConfig{Float64}()

"""
    BuddyCheckResult{T<:AbstractFloat}

Result of buddy check for a single observation.

# Fields
- `passed::Bool`: Whether observation passed buddy check
- `failure_reason::String`: Reason for failure (if failed)
- `num_buddies::Int`: Number of buddy observations found
- `buddy_mean::T`: Mean value of buddy observations
- `buddy_std::T`: Standard deviation of buddy observations
- `normalized_departure::T`: Normalized departure from buddy mean
- `qc_flag::Int`: QC flag (0=good, 1=suspect, 2=bad)
"""
struct BuddyCheckResult{T<:AbstractFloat}
    passed::Bool
    failure_reason::String
    num_buddies::Int
    buddy_mean::T
    buddy_std::T
    normalized_departure::T
    qc_flag::Int
end

"""
    perform_buddy_check(obs::AbstractObservation{T},
                       all_observations::Vector{<:AbstractObservation{T}},
                       config::BuddyCheckConfig{T}) where T

Perform buddy check on observation against nearby observations.

# Arguments
- `obs`: Observation to check
- `all_observations`: All available observations for buddy comparison
- `config`: Buddy check configuration

# Returns
- `BuddyCheckResult{T}`: Detailed check results

# Examples
```julia
config = BuddyCheckConfig{Float64}()
result = perform_buddy_check(obs, all_obs, config)
if result.passed
    println("Observation passed buddy check")
    println("  Buddies found: ", result.num_buddies)
    println("  Buddy mean: ", result.buddy_mean)
else
    println("Rejected: ", result.failure_reason)
end
```
"""
function perform_buddy_check(obs::AbstractObservation{T},
                             all_observations::Vector{<:AbstractObservation{T}},
                             config::BuddyCheckConfig{T}) where T

    # Find buddy observations
    buddies = find_buddy_observations(obs, all_observations, config)

    # Check if sufficient buddies found
    if length(buddies) < config.min_buddies
        # Not enough buddies - pass by default but flag
        return BuddyCheckResult{T}(
            true,  # Pass (insufficient data to reject)
            "insufficient_buddies",
            length(buddies),
            obs.observed_value,  # Use own value as buddy mean
            zero(T),
            zero(T),
            0  # Good (but noted)
        )
    end

    # Compute buddy statistics
    buddy_mean, buddy_std = compute_buddy_statistic(obs, buddies, config)

    # Compute normalized departure
    departure = obs.observed_value - buddy_mean

    # Account for observation error and buddy uncertainty
    buddy_uncertainty = sqrt(
        obs.observation_error.total_error^2 +
        (buddy_std^2 / length(buddies))
    )

    normalized_departure = abs(departure) / buddy_uncertainty

    # Apply buddy check threshold
    passed = normalized_departure <= config.rejection_threshold

    qc_flag = if passed
        0  # Good
    elseif normalized_departure <= config.rejection_threshold * 1.5
        1  # Suspect
    else
        2  # Bad
    end

    failure_reason = passed ? "" : "buddy_check_failed"

    return BuddyCheckResult{T}(
        passed,
        failure_reason,
        length(buddies),
        buddy_mean,
        buddy_std,
        normalized_departure,
        qc_flag
    )
end

"""
    find_buddy_observations(obs::AbstractObservation{T},
                           all_observations::Vector{<:AbstractObservation{T}},
                           config::BuddyCheckConfig{T}) where T

Find suitable buddy observations for comparison.
"""
function find_buddy_observations(obs::AbstractObservation{T},
                                all_observations::Vector{<:AbstractObservation{T}},
                                config::BuddyCheckConfig{T}) where T

    buddies = AbstractObservation{T}[]

    # Get search radius for this variable type
    search_radius = get(config.variable_specific_radius, obs.variable, config.horizontal_radius)

    for other_obs in all_observations
        # Skip self
        if other_obs === obs
            continue
        end

        # Must be same variable type
        if other_obs.variable != obs.variable
            continue
        end

        # Check horizontal distance
        distance = geographic_distance(obs.location, other_obs.location)
        if distance > search_radius
            continue
        end

        # Check temporal difference
        time_diff = abs(Dates.value(other_obs.observation_time - obs.observation_time)) / 1000.0  # seconds
        if time_diff > config.temporal_window
            continue
        end

        # Check vertical separation (for upper-air observations)
        if needs_vertical_check(obs)
            vertical_sep = abs(get_pressure(other_obs) - get_pressure(obs))
            if vertical_sep > config.vertical_radius
                continue
            end
        end

        # Check elevation difference (for surface observations)
        if config.elevation_correction && is_surface_observation(obs)
            elev_diff = abs(other_obs.location.altitude - obs.location.altitude)
            if elev_diff > 500.0  # 500 meters max elevation difference
                continue
            end
        end

        # Observation passed all criteria - add as buddy
        push!(buddies, other_obs)
    end

    return buddies
end

"""
    compute_buddy_statistic(obs::AbstractObservation{T},
                           buddies::Vector{<:AbstractObservation{T}},
                           config::BuddyCheckConfig{T}) where T

Compute buddy mean and standard deviation with optional weighting.
"""
function compute_buddy_statistic(obs::AbstractObservation{T},
                                buddies::Vector{<:AbstractObservation{T}},
                                config::BuddyCheckConfig{T}) where T

    if isempty(buddies)
        return obs.observed_value, zero(T)
    end

    buddy_values = [b.observed_value for b in buddies]

    if !config.distance_weighting
        # Simple arithmetic mean
        buddy_mean = mean(buddy_values)
        buddy_std = std(buddy_values; corrected=true)
        return buddy_mean, buddy_std
    end

    # Distance-weighted mean
    distances = [geographic_distance(obs.location, b.location) for b in buddies]
    buddy_mean, buddy_std = distance_weighted_mean(buddy_values, distances, config)

    # Apply elevation correction if needed
    if config.elevation_correction && is_surface_observation(obs)
        buddy_mean = apply_elevation_correction(
            obs, buddy_mean, buddies, config
        )
    end

    return buddy_mean, buddy_std
end

"""
    distance_weighted_mean(values::Vector{T}, distances::Vector{T}, config::BuddyCheckConfig{T}) where T

Compute distance-weighted mean and standard deviation.

Uses exponential weighting: w(d) = exp(-d/L) where L is length scale.
"""
function distance_weighted_mean(values::Vector{T}, distances::Vector{T},
                               config::BuddyCheckConfig{T}) where T

    if isempty(values)
        return zero(T), zero(T)
    end

    # Length scale for exponential weighting (1/3 of search radius)
    length_scale = config.horizontal_radius / 3.0

    # Compute weights
    weights = exp.(-distances ./ length_scale)
    weights ./= sum(weights)  # Normalize

    # Weighted mean
    weighted_mean = sum(weights .* values)

    # Weighted standard deviation
    weighted_var = sum(weights .* (values .- weighted_mean).^2)
    weighted_std = sqrt(weighted_var)

    return weighted_mean, weighted_std
end

"""
    apply_elevation_correction(obs::AbstractObservation{T},
                              buddy_mean::T,
                              buddies::Vector{<:AbstractObservation{T}},
                              config::BuddyCheckConfig{T}) where T

Apply elevation correction to buddy mean for surface observations.

For temperature and pressure, adjusts for elevation differences using
standard atmospheric lapse rate.
"""
function apply_elevation_correction(obs::AbstractObservation{T},
                                   buddy_mean::T,
                                   buddies::Vector{<:AbstractObservation{T}},
                                   config::BuddyCheckConfig{T}) where T

    if !is_surface_observation(obs) || isempty(buddies)
        return buddy_mean
    end

    # Compute mean buddy elevation
    mean_buddy_elev = mean([b.location.altitude for b in buddies])
    elev_diff = obs.location.altitude - mean_buddy_elev

    # Apply variable-specific correction
    corrected_mean = buddy_mean

    if obs.variable == TEMPERATURE
        # Standard lapse rate: -6.5 K/km
        lapse_rate = T(-6.5)  # K/km
        corrected_mean = buddy_mean + (elev_diff / 1000.0) * lapse_rate

    elseif obs.variable == SURFACE_PRESSURE
        # Hydrostatic correction: dp/dz = -ρg ≈ -p/(RT) * g
        # Approximate: Δp/p ≈ -Δz/H where H ≈ 7400 m
        scale_height = T(7400.0)  # meters
        pressure_ratio = exp(-elev_diff / scale_height)
        corrected_mean = buddy_mean * pressure_ratio

    elseif obs.variable == DEW_POINT
        # Dew point typically decreases with altitude (similar to temperature)
        # Use similar lapse rate but slightly reduced
        lapse_rate = T(-5.0)  # K/km
        corrected_mean = buddy_mean + (elev_diff / 1000.0) * lapse_rate
    end

    return corrected_mean
end

# Helper functions

"""
    geographic_distance(loc1::ObservationLocation{T}, loc2::ObservationLocation{T}) where T

Compute great circle distance between two locations using haversine formula.
Returns distance in meters.
"""
function geographic_distance(loc1::ObservationLocation{T}, loc2::ObservationLocation{T}) where T
    # Earth radius in meters
    R = T(6371000)

    # Convert to radians
    lat1_rad = deg2rad(loc1.latitude)
    lat2_rad = deg2rad(loc2.latitude)
    lon1_rad = deg2rad(loc1.longitude)
    lon2_rad = deg2rad(loc2.longitude)

    # Haversine formula
    dlat = lat2_rad - lat1_rad
    dlon = lon2_rad - lon1_rad

    a = sin(dlat/2)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(dlon/2)^2
    c = 2 * atan(sqrt(a), sqrt(1-a))

    return R * c
end

"""
    needs_vertical_check(obs::AbstractObservation) -> Bool

Determine if observation type requires vertical buddy check.
"""
function needs_vertical_check(obs::AbstractObservation)
    return obs isa RadiosondeObservation ||
           obs isa AircraftObservation ||
           obs isa GPSROObservation
end

"""
    is_surface_observation(obs::AbstractObservation) -> Bool

Determine if observation is a surface observation.
"""
function is_surface_observation(obs::AbstractObservation)
    return obs isa SurfaceObservation
end

"""
    get_pressure(obs::AbstractObservation{T}) where T -> T

Extract pressure level from observation.
"""
function get_pressure(obs::RadiosondeObservation{T}) where T
    return obs.pressure_level
end

function get_pressure(obs::AircraftObservation{T}) where T
    return obs.location.pressure
end

function get_pressure(obs::GPSROObservation{T}) where T
    return obs.location.pressure
end

function get_pressure(obs::AbstractObservation{T}) where T
    return obs.location.pressure
end

"""
    platform_quality_weight(obs::AbstractObservation{T}) where T -> T

Compute platform-dependent quality weight for buddy check.

Higher quality platforms get more weight in buddy calculations.
"""
function platform_quality_weight(obs::AbstractObservation{T}) where T
    if obs isa RadiosondeObservation
        return T(1.0)  # Highest quality
    elseif obs isa SurfaceObservation
        if obs.station_type == "SYNOP"
            return T(0.9)
        else
            return T(0.8)
        end
    elseif obs isa AircraftObservation
        if obs.data_source == "AMDAR"
            return T(0.85)
        else
            return T(0.75)
        end
    elseif obs isa RadianceObservation
        return T(0.7)  # Lower weight due to representation error
    else
        return T(0.5)  # Unknown platform
    end
end

end # module BuddyCheck
